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Abstract. We report a first principles study of the coupled evolution of energetic 
' ions, background majority ions, electrons and electromagnetic fields in magnetised 



plasma during the linear phase of the lower hybrid drift instability. A particle-in- 
cell code, with one spatial and three velocity space co-ordinates, is used to analyse the 
O ' evolving distribution of a drifting ring-beam population of energetic protons in physical 

c/3 i space and gyrophase angle. This analysis is carried out for plasma parameters that 

approximate to edge conditions in large tokamaks, in a scenario that is motivated 
by observations of ion cyclotron emission and may be relevant to alpha channelling. 
Resonant energy transfer occurs at the two gyrophase angles at which the instantaneous 
speed of an energetic proton on its cyclotron orbit precisely matches the phase 
■ velocity of the lower hybrid wave along the simulation domain. Electron space- 

CTN , charge oscillations determine the wavelength of the propagating lower hybrid wave, 

and thereby govern the spatial distribution of gyrobunching of the energetic protons 
. that drive the instability. 



1. Introduction 

^ . Understanding the nature of the wave-particle interaction that drives the lower 
hybrid drift instability (LHDI) p3-H] is of fundamental interest in plasma physics. 
The LHDI is potentially operative wherever drifting populations of ions (whether 
suprathermal minorities, or arising from equilibrium gradients) occur in magnetised 
plasmas. There is an extensive literature documenting its many roles in both space [HH7] 
and, increasingly, laboratory |8, 9j plasma physics. Predominantly electrostatic waves, 
propagating quasi-perpendicular to the magnetic field, are preferentially excited by 
the LHDI. Their frequency is characterised by the lower hybrid frequency u>lh = 
(f2 ce f2 ci /(l + riceflci/uipi)) 1 / 2 , which exceeds the cyclotron frequency of the ions and is 
below that of the electrons. It is clear from the foregoing that the LHDI possesses 
features that invite deeper investigation. The instability arises from the beam-like (i.e., 
directed in velocity space) character of the fast ion population that drives it, while 
the waves excited are conditioned by the cyclotronic (i.e., rotational in velocity space) 
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natural frequency Vt ci of the background thermal ions. The waves excited are conditioned 
by the electron- ion mass ratio through Q ci /Q ce , implying that electron inertia plays a 
role. The waves excited are also conditioned by the density of charge carriers through 
the ion plasma frequency u P i, implying that charge separation plays a central role. 

In the present paper we study the physics of the excitation mechanism of the 
LHDI, using a particle-in-cell (PIC) code [TOl - TH] which provides a fully kinetic 
description of the interaction between minority (~ 1%) energetic protons, background 
thermal deuterons, electrons, the self-consistent electromagnetic fields, and the imposed 
background magnetic field. In particular we focus on the dynamics of the ions as they 
drive the LHDI, which is found to exhibit coherent bunching that is highly structured 
in both velocity space and real space. To anticipate some key results, the kinetic 
treatment of both electrons and ions enables us to capture two essential features of 
the physics that are not addressed (because they are averaged over) in higher order 
reduced models. First, the code directly represents the excited waves as propagating 
electron charge density oscillations with well defined frequency and wavelength. Second, 
the code directly represents ion motion and can therefore resolve the angular sector of an 
ion cyclotron orbit during which the ion is transiently in Doppler-shifted phase resonance 
with the wave that is being excited. 

The scenario for the LHDI considered here is motivated by observations of radiative 
collective instability (ion cyclotron emission, ICE) of minority fusion-born and beam- 
injected energetic ions in the outer mid-plane edge plasma of the JET [ToTflT] and TFTR 
[T8H20] tokamaks. Analysis of the spatial drift orbits of the ions responsible for exciting 
the fast Alfven waves observed showed [T5l[T9] that their velocity space distribution 
can be approximated by the analytical model f(v±,v\\) = l/(2irv r )5(v± — v r )5(v\\ — it) 
where v r and u correspond to a pitch angle just inside the trapped-passing boundary 
for an ion near its birth energy. In Refs. [2Tjl22j it was shown that these distributions 
can also be unstable against the LHDI. Furthermore the lower hybrid waves excited 
were found to undergo Landau damping on resonant electrons asymmetrically with 
respect to the distribution of parallel velocities, resulting in the spontaneous creation 
of a current supported by suprathermal electrons J2TJ|22]. This is a key building block 
of alpha channelling [23|l24] scenarios for fusion plasmas. The principle of lower hybrid 
current drive by Landau damping of externally applied waves was identified by Fisch 
and colleagues [25j|26] for fusion plasmas, and an application to spacecraft and rocket- 
borne measurements of electron distributions in Earth's auroral zone was noted in 
Ref. [27J. The JET and TFTR fusion plasma results motivate the present choice of: 
plasma parameters, approximating to tokamak edge conditions, where computationally 
affordable; and energetic ion velocity distribution, relevant to ICE observations [T5TJ20] 
and to recent simulations of LHDI alpha channelling [2TJ[22]. While these represent 
only one possible realisation of an LHDI scenario, among many, we believe that the 
underlying physics of the phase resonant excitation described in the present paper is 
generic to a significant extent. 
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2. Fundamentals of the simulation 



Figure 1 depicts the unperturbed spatiotemporal orbit, during one cyclotron period, of 
a single energetic proton in a field-aligned coordinate system, and a coordinate system 
which is tilted with respect to the magnetic field. The initial and final projection of the 
equilibrium cyclotron orbit on the plane perpendicular to B is also shown for guidance. 
This figure shows how the energetic ion distribution f(v±,v\\) ~ S(v± — v r )6(v\\ — u) 
governs the movement of each ion along the simulation domain x, and the gyrophase 
a = arctan(f_|_ jl /v_|_ 2 ). There are two directions of spatial variation in two co-ordinate 
systems shown in the figure. Only one of these is the single spatial direction x of the 
PIC code; the other three are for schematic purposes. The x-projected motion is shown 
as a function of gyrophase a in Fig. 2, which implies that the particle velocity v x = dx/dt 
along the simulation domain varies continuously, so that phase resonance between the 
particle v x and any wave that can propagate with phase velocity u/k x will only occur 
transiently, along a particular segment of its cyclotron orbit. It is in this segment that 
the perturbed (x, a) phase orbit of the ion resonantly exciting the wave will deviate 
most strongly from the unperturbed orbit of Fig.2. 



Figure 1. Unperturbed motion through time (indicated by shading (colour online)) 
of an energetic proton, plotted with respect to both magnetic field-based (upright solid 
cube) and Cartesian (tilted dashed cube) spatial co-ordinates. The plot follows one 
cyclotron period and incorporates drift antiparallel to B, which is inclined at = 84° 
to the one-dimensional spatial simulation domain x. Here r_\_,x and r± t 2 are orthogonal 
co-ordinates perpendicular to B. The gyrophase a is zero when the proton passes 
through the plane defined by x and B. Projections of the cyclotron motion on the 
plane perpendicular to B are shown at t = and one cyclotron period later. The 
combination of cyclotron motion and drift parallel to B results in the drift shown in 




Fig.2. 
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Figure 2. Projected motion through time (indicated by colour) of an energetic proton 
along the simulation domain direction x as a function of gyrophase a. This motion is 
implicit in the definition of the unperturbed velocity distribution in terms of v± and 
v\\ , and follows from Fig.l. The parallel velocity component present in this geometry 
causes a drift in x of ut cp cos9 per gyro-orbit, where r cp — 27r/f2 cp is the proton 
gyroperiod. 

It is known from Figs. 3 and 4 of Ref. |21j that the dominant collectively excited 
modes propagate in the negative x-direction with u/k x ~ 1.5 x lCr 7 m/s < v x = 
ucos6 — v r sind. Modes with slightly lower amplitude are excited propagating in the 
opposite direction with similar phase speed. These phase velocities are located near 
the maxima of transient speeds, in either direction along the simulation domain, of the 
energetic ions that drive the instability. These occur (see Figs. 1 and 2) at gyrophase 
a ~ it/2 and ~ 3tt/2; see also panel (a) of Fig. 2 of Ref. [22J. 

3. Simulation results 

We now turn to the resonant wave-particle interaction as manifested by the distribution 
of energetic protons along the physical co-ordinate x and with respect to gyrophase 
angle a. Figure 3 plots (left) the whereabouts of particles in (x, a) phase space at 
t = \9.8tlh, towards the end of the linear phase of the LHDI in this scenario, and 
(right) the corresponding distribution of particle energies. Within the spatial limits of 
the simulation box, there are 21 gyrophase bunched structures which span the whole 
gyrophase range: 10.5 of these are visible in Fig 3, because data from only half of the 
simulation domain is shown (Figs 3 to 6 all show data from only half of the spatial 
domain for reasons of pictorial resolution). The simulation domain length L ~ 6.3 
proton Larmor radii. This is sufficiently large that protons undergo ~ 7 gyro-orbits 
per boundary crossing, which corresponds to approximately one boundary crossing per 
proton during the entire simulation. Figure 3 demonstrates that the rapidly varying 
behaviour of the distribution in the negative a domain of phase space is a continuation 
of the slower variation in the positive domain. This new result is not apparent from 
previous analysis of the LHDI in this regime, for example Fig. 5 of Ref. [21] and Fig.4 of 
Ref. [22]. 
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Figure 3. PIC code derived proton population data from a snapshot at t = 19.8t£h 
towards the end of the linear stage of the instability, plotted as a function of gyrophase 
(ordinate axis) and position (abscissa) . Position axes are in units of box length L and 
gyrophase axes are in radians. Panel (a): Proton probability density on a logio scale 
indicated by shading (colour online). Panel (b): Proton energy in units of e p (t = 0) 
indicated by shading (colour online). Both panels show 10.5 nearly-identical s-shaped 
features in the upper half plane at regular intervals in the spatial direction, as well as 
smaller scale repeating structures. 

There are 10.5 gyrophase bunched structures per half-box because the simulation 
box accommodates 21 wavelengths of the wave that is found to be spontaneously excited. 
This is demonstrated in Fig. 4, which plots the spatial distribution of electrons at the 
snapshot in time shown in Fig. 3, and also confirms that the predominantly electrostatic 
excited wave is primarily a charge density oscillation supported by the electrons. 




Figure 4. PIC code derived electron population data from a snapshot at t = 19.8t£# 
towards the end of the linear stage of the instability. Electron probability density on a 
logio scale is indicated by shading (colour online) as a function of gyrophase (ordinate 
axis) and position (abscissa). Position axes are in units of box length L and gyrophase 
axes are in radians. The vertical stripes indicate that electron bunching occurs as a 
function of position and not gyrophase: electrons participate in wave action by charge 
density oscillation. 
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Figure 3 (left) captures wave excitation by the energetic ions as it occurs. The 
strong spatially dispersive distortion at a = it /2 corresponds (see Fig.l) to the maximum 
negative particle velocity (see also Fig. 2 (a) of Ref. [22]), which is located between the 
points in gyrophase where the resonant drive of the dominant backward-propagating 
waves takes place. The speed of the dominant wave is slightly smaller than the maximum 
speed of the protons. Consequently there are two positions of resonance in gyrophase, 
either side of the maximum negative velocity. These resonant perturbations apparently 
condition the long-wavelength (in gyrophase a) structure across the positive-a domain 
of phase space. At a = — n/2, which corresponds to the maximum positive particle 
velocity, the consequences of resonance are apparent in the moire pattern (which is 
more easily visible in Fig.3 (right)). While both these velocity resonances give rise to 
significant (but differing by an order of magnitude) concentrations of electromagnetic 
field energy in (u, fc)-space, as shown in Fig.3 of Ref. [21J, it follows from Fig.3 that the 
phase space consequences are different in the two cases. The number of moire s-shaped 
features in the lower half plane of Fig.3 exceeds that in the upper half plane. This 
reflects the difference in k x for the waves excited at the two resonances, see again Fig. 
3 of Ref. [21]. 

Figure 3 (right) assists interpretation of the physics in terms of the simplified 
analytical model for the dynamics of the particle-wave interaction that was introduced 
in Eqs.(12) to (16) of Ref. [22]- The essence of this model is that the impact of wave- 
particle interaction, for a proton with initial gyrophase ciq at t = and at arbitrary (x, t), 
scales with the magnitude of the electric field which that proton experienced when it 
was most closely in phase resonance with the wave. This resonance will have occurred 
at time £r(oio) such that v x (tji) = u/k, where oj ~ Gujlh and k = —2u pe /c are the 
angular frequency and wavenumber of the dominant excited mode in the simulation, 
obtained from the analysis of the Fourier transform of the electromagnetic field in 
Figs. 3 of Refs. [21J. Given an energetic proton moving on an unperturbed orbit drawn 
from the drifting beam distribution specified previously, expressions for t^ and position 
xn(aQ,x,t) follow from Eqs.(14) and (15) of Ref. [22]- The key point is that the field 
amplitude, ^(ao,x, t) experienced at the point of resonance (xji,t R ), which is given 
by Eq.(16) of Ref. [22J, is considered a property of the proton at all subsequent (x,t). 
Figure 5 (left) plots S^cto^x, t) for a population of energetic protons that is initially 
uniformly distributed in space x and gyrophase angle a at an arbitrary time t = 2ti /Q cp . 
Figure 5 (right) is derived from the analytical model using the three waves identified 
from the Fourier transform of the excited electromagnetic field that is generated by 
the PIC simulation during its linear phase. To construct it we repeat the operation 
used for Fig. 5 (left) while adding: a second wave whose (u, k) are those identified 
as the harmonic of the dominant backward travelling wave in the Fourier transform 
{ojh = !2uLH,kh = —4:U pe /c), and whose amplitude is smaller by a factor ten; and a 
third wave whose (u, k) are those identified for the forward travelling excited wave in 
the Fourier transform (cjf = 6.6ulh, kf = l-3w pe /c), and whose amplitude is smaller by 
a factor four. 
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Figure 5 (right) appears to capture the key features of Fig. 3 (right). This analytical 
model is based on unperturbed proton gyro-orbits, and Fig. 5 shows that these suffice 
to recreate gyrobunching effects seen in the PIC simulations provided perturbations 
are small in the sense that: the resonant gyrophase angles are sharply defined; and 
the spatial deviation of the resonant particles from their unperturbed orbits is small 
compared to the wavelength. 

This analytical model is based on unperturbed proton gyro-orbits and is able to 
recreate gyrobunching effects seen in the PIC simulation. This is the case because 
perturbations are small, in the sense that the resonant gyrophase angles are sharply 
defined and the spatial deviation (as distinct from gyrophase deviation) of the resonant 
particles from their unperturbed orbits is also small. 




x[L] x[L] 

(a) (b) 

Figure 5. Normalized electric field amplitude S indicated by shading (colour online) 
experienced by protons at their most recent point of resonance, as a function of proton 
gyrophase and position at a snapshot in time (see text). Panel (a): Normalized wave 
amplitude seen by protons at resonance with the dominant wave. Panel (b): As panel 
(a) for the sum of dominant wave of unit amplitude, its second harmonic of one tenth 
amplitude, and a wave travelling in the opposite direction of one quarter amplitude. 
Data in panel (a) shows the effect of resonance with the dominant exited wave as 
visible from the large scale structure in Fig. 3a and b, while panel (b) shows finer 
scale structures, as in Fig. 3b, highlighting the effects of resonances with less powerful 
excited waves. 

4. Conclusions 

Analysis of the distribution in physical space and gyrophase angle of a drifting ring- 
beam population of energetic ions, as they drive the linear phase of the lower hybrid 
drift instability in a 1D3V PIC simulation, yields insights into the underlying plasma 
physics. The spatial coherence and spatial variation of the gyrophase bunched structures 
are shown to reflect a combination of cyclotron-type and beam-type dynamical and 
resonant effects, conditioned by the interplay of ion and electron dynamics. These new 
results shed fresh light on the hybrid character of the lower hybrid drift instability. 
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Resonant energy transfer occurs at the narrowly defined gyrophases at which the 
instantaneous speed of an energetic proton, from the drifting ring-beam population, 
on its cyclotron orbit precisely matches the phase velocity of the lower hybrid wave 
along the simulation domain. As in any resonant wave-particle instability, coupling to 
a wave that can propagate is necessary. In the present case, it is electron space-charge 
oscillation (Fig. 4) that determines the wavelength of the propagating lower hybrid wave, 
and thereby governs the spatial distribution of gyrobunching of the energetic protons 
(Figs. 3 and 5) that drive the instability. The gyrophase structure of phase space 
bunching at any given position is governed by the proton trajectories on their cyclotron 
orbits, which are oblique to the direction of propagation of waves along the simulation 
domain, in conjunction with restrictions on the spatial separation of gyrobunches due 
to collective electron oscillation. 

Many of the key features of the wave-particle interaction that we observe in the 
simulations can be understood retrospectively in terms of a simplified analytical model. 
However this model is wholly reliant on PIC code outputs for its specification, including 
the number of waves, their relative amplitudes, and their frequencies and wavenumbers. 

This analysis was carried out for plasma parameters that approximate, where 
computationally affordable, edge plasma conditions in a large tokamak. While the 
physics is generic, we note that this scenario is motivated by observations of ion cyclotron 
emission from JET [15J and TFTR [19] , and may be relevant to alpha channelling [211122] . 
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